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Abstract 

Aeroassisted planetary entry uses atmospheric drag to decelerate spacecraft from super-orbital to orbital or sub- 
orbital velocities. Numerical simulation of flow fields surrounding these spacecraft during hypersonic atmospheric 
entry is required to define aerothermal loads. The severe compression in the shock layer in front of the vehicle and 
subsequent, rapid expansion into the wake are characterized by high temperature, thermo-chemical nonequilibrium 
processes. Implicit algorithms required for efficient, stable computation of the governing equations involving 
disparate time scales of convection, diffusion, chemical reactions, and thermal relaxation are discussed. Robust 
point-implicit strategies are utilized in the initialization phase; less robust but more efficient line-implicit strategies 
are applied in the endgame. Applications to ballutes (balloon-like decelerators) in the atmospheres of Venus, Mars, 
Titan, Saturn, and Neptune and a Mars Sample Return Orbiter (MSRO) are featured. Examples are discussed where 
time-accurate simulation is required to achieve a steady-state solution. 
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total energy per unit mass 
inviscid flux vector 

viscous flux vector 
total enthalpy 

cell index or reference length 
Mach number 
mixture molecular weight 
molecular weight of species i 
inviscid flux part of SI L , Eq. (2) 
viscous flux part of Eq. (3) 
pressure, N/m I 2 
dependent variable vector 
wall heat transfer rate, W/cm 2 
residual vector at cell L 
universal gas constant 

, Reynolds number 

source term part of 3$ L 
translational-rotational temperature, K 
vibrational-electronic temperature, K 
velocity component a direction, (X=x, y, or z 
free stream velocity 
Cartesian coordinates 


* Senior Research Engineer, Aerothermodynamics Branch, Associate Fellow AIAA 

Copyright © 2001 by the American Institute of Aeronautics and Astronautics, Inc. No copyright is asserted in the United States 
under Title 17, U.S. Code. The U.S. Government has a royalty-free license to exercise all rights under the copyright claimed 
herein for Governmental Purposes. All other rights are reserved by the copyright owner. 


I 

American Institute of Aeronautics and Astronautics 



^ ro 


mole fraction of species i 
eigenvalue limiting factor, 0<£<0.3 
eigenvalue 
|i viscosity 

p density 

a cell wall area 

Q cell volume 

Subscripts 

i species index 

m cell wall index 

w wall surface conditions 

oo free stream conditions 


from the point-implicit formulation in the center 
diagonal but provide enhanced convergence. 

This paper reviews these implicit algorithms 
as applied in the Langley Aerothermodynamic 
Upwind Relaxation Algorithm (LAURA) 2 to a 
variety of challenging (from perspective of physical 
and geometric complexity) problems in aeroassist 
simulation. It defines performance benchmarks for a 
model problem so that other algorithms may be 
compared and more efficient simulation tools may 
evolve. Finally, it notes examples of a commonly 
applied, local-time stepping (constant Courant 
number) procedure which precludes attainment of a 
steady state solution in some cases where a steady 
state exists. 


Introduction 

Aeroassisted planetary entry uses 
atmospheric drag to decelerate spacecraft from super- 
orbital to orbital or sub-orbital velocities. The 
vehicles encounter flow domains from free-molecular 
to continuum 1 . The severe compression in the shock 
layer in front of the vehicle and subsequent, rapid 
expansion into the wake are characterized by high 
temperature, thermo-chemical nonequilibrium 
processes. The chemical state of atmospheric 
constituents crossing the bow shock varies from 
frozen to nonequilibrium (ionization, dissociation) to 
equilibrium. The thermal state of internal energy 
modes (electronic, vibrational, rotational) varies 
similarly. Convective heating rates are strongly 
influenced by the catalytic properties of the surface 
and diffusion rate of atoms to the wall. 

Numerical simulations are the primary 
source for defining environments and vehicle 
response to environments required in the design of 
planetary entry vehicles. The simulations are verified 
when possible with ground-based and flight test data 
that cover some subset of the entry corridor 
parameter space. Implicit algorithms are required for 
efficient, stable computation of the governing 
equations involving the disparate time scales of 
convection, diffusion, chemical reactions, and 
thermal relaxation. Point- and line-implicit 
treatments provide opportunities for single level 
storage relaxation at the expense of time accuracy. 
The essence of a point-implicit method is to treat all 
contributions to the residual from the cell center at 
the advanced time level — upwind formulations 
enable this structure by introducing cell center 
dependence from the convective terms. Implicit lines 
running across the viscous layer, solved by block tri- 
diagonal relaxation, utilize the same Jacobian blocks 


Algorithm 

Point-Implicit Relaxation: 

The Langley Aerothermodynamic Upwind 
Relaxation Algorithm (LAURA) employs two 
relaxation modes: point-implicit and line-implicit. 

The point- implicit mode 2,3 implicitly solves only 
contributions of dependent variables at cell L, q, , to 
the residual at cell L, 7, ; contributions from 
dependent variables at neighbor cells are treated 
explicitly. The contributions of inviscid, viscous, and 
source terms are included in the point- implicit 
matrix, as follows: 

= S[Q, + CjJij + c 2 jV, (1) 

where S L is the Jacobian of the chemical and thermal 
source term components of the residual, is the 
Jacobian of the inviscid terms, and jV l is the Jacobian 
of the viscous terms - all with respect to q, . Scalar 
constants d>1.5 and c 2 >0.5 are relaxation factors. 
Typical values are c ( =3 and c 2 =l. Relaxation factors 
on the source term Jacobian have never been 
required. 

Within the context of flux difference 
splitting (FDS) and Roe's averaging 4 the inviscid 
Jacobian is defined: 

•*.. = ^XI A J o " (2) 

where m is summed over all faces of cell L, A m is the 
Roe’s averaged flux Jacobian at face m, and |A ml ' S 

defined using the absolute values of the eigenvalues 
of A m . The elements of A m are required in the 
definition of the FDS formulation of the inviscid 
flux; consequently, it is economical to reuse these 
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quantities in the formulation of the Jacobian. The 
absolute value of the eigenvalues may be limited 
from below using a variation of Harten’s entropy 
fix 5 ’ 3 . This limiter serves to enhance stability of the 
relaxation algorithm and prevent formation of 
expansion shocks across sonic lines or carbuncles in 
the stagnation region. The limiter is defined such 
that it is neither too large - adversely effecting 
computed surface heating rates and shear forces, nor 
too small - adversely effecting stability. In the 
direction tangent to the boundary layer it is limited by 
normal to the boundary layer it is limited by 

£ ■ maxjVu 2 + v 2 + w" , 0.0001 • a . 

The viscous Jacobian is defined: 

A = X B m c m (3) 


where B m is the Jacobian of viscous terms across the 
m* face of cell L. It implicitly treats only the divided 
difference across the cell face. Transport properties at 
cell face m are linear averages of adjacent cell center 
values. Velocities in the shear work term are the Roe 
averaged quantities already computed for the inviscid 
flux function. Eigenvalues of B m are positive or zero. 

The change in dependent variables after a 
single, point-implicit relaxation step is computed as 
follows: 


accommodates a very simple, self-contained 
initialization procedure for a generic family of 
spacecraft geometries that generates grids in a thin 
layer normal to the surface and allows the grid to 
grow out as the solution evolves. It easily handles 
multi-block, asynchronous relaxation. The only 
numerical “knobs” to control the relaxation process 
are the magnitudes of the relaxation factors c 5 and c 2 , 
the frequency of Jacobian updates, and the frequency 
of grid re-alignments. 


Line-Implicit Relaxation: 

The point-implicit option is not efficient for 
relaxing the viscous layers. Some efficiency can be 
recovered by using multitasking with majority of 
cycles in the boundary layer or by using an over- 
relaxation factor 1>c 2 >0.5. However, a better 
approach is to introduce implicit dependence of the 
solution across the entire shock layer. 

The line-implicit relaxation algorithm is a 
simple extension of the point-implicit algorithm. 
However, approximations to (X the Jacobian of 
r, with respect to q l _ l , and C L , the Jacobian of 
r L with respect to q, +1 , are introduced to preserve 


symmetry of form with that also improves 
convergence. They are defined: 


^ = 



Aj.-w 2 ) 


— c 2 B, _ 1/2 


(5a) 


Q. 

+ a, 


Aq, = r L 


(4) 


where I - the identity matrix, Q L - the volume, and 
At, - the time step at cell L (usually derived from a 

constant Courant number specification) are 
associated with the time derivative term of the 
conservation equations. LAURA generally uses a 
Courant number of one million (a philosophical 
constraint to admit only algorithms that will solve a 
steady equation system without dependence on time 
step); consequently, this term may be neglected 
relative to «® L . In cases with crudely defined initial 
conditions or with difficult transients stability is 
maintained by temporarily increasing the relaxation 
factors ci and c 2 . Values of 10 are usually sufficient; 
values up to 100 have been required on some highly 
energetic, multi-species aeroassist simulations. 

The point-implicit option in LAURA is 
generally robust - though certainly not a black box. 

It works welt in conjunction with a one-dimensional 
grid adaptation algorithm in LAURA that makes 
large-scale adjustments of the inflow boundary to 
align with evolving bow shock location. It 


^ = 


^ (Ja ,+ 1 / 2 1 a, +1/2 ) 

— C 2® L+l / 2 


(5b) 


It is assumed that sequential indices L define a 
continuous series of cells spanning some portion of 
the computational domain. The subscripts L ± 1 / 2 
refer to Roe averaged conditions on face m between 
cell centers L and L ± 1 , respectively. A more exact 
linearization of the residual would lead to use of 
Aj ±1 in place of A, ±l/2 in the second occurrence of 
Eqs. 5; however, the present formulation provides 
exact cancellation of convective, downwind implicit 
influences and exhibits more robust convergence 
characteristics, particularly in the vicinity of strong 
shock waves at large Courant numbers. 

In the present work, the source term requires 
no off-diagonal, implicit contribution. The line- 
implicit relaxation in Eq. (6) requires approximately 
2.8 times more CPU time per iteration as the point- 
implicit algorithm in the applications described here 
but convergence rate more than makes up for this 
overhead. It requires additional storage for the off- 
diagonal blocks. 


3 

American Institute of Aeronautics and Astronautics 


^1+1 

0 

0 

0 


e, o o o 

®i.+i +i ® ® 

0 


0 ^»,+ N-l ®I.+ N-I ^l.+ N-l 

0 0 ^l.+ N ®I. + N 


1 

< 

1 


r. 

Aq,.+i 


?L+I 

Aq,. tN -i 


r i.+ N-l 

. Aq,. + N . 


r i.+ N 


( 6 ) 


Implicit Boundary Conditions: 

Implicit treatment of boundary conditions is 
required to obtain good convergence with line 
relaxation. Boundary conditions in the present 
formulation are implemented with ghost cells. If the 
off-diagonal Jacobian at the ghost cell for L= 1 is flo 
then is augmented by (LqWq A( U) x where 

n dq„ = ^dq, (7) 

is an expression of the dependence of the ghost cell 
on the boundary cell. In like manner, <$ N is 
augmented by Cn+iWn+i^n where 
W N+l dq N + l = W N dq N . The only non-trivial 
boundary condition considered here (see Appendix 
A) is at solid walls that account for finite catalysis 
and radiative equilibrium wall temperatures. 

Time Accuracy: 

The LAURA algorithm, in its default mode, 
uses single level storage and relaxes the steady form 
of the governing equations with a Courant number of 
10 6 . However, some applications considered here 
require time-accurate relaxation. Time accuracy 
(first-order in At ) can be recovered by using constant 
time step advancement at all cells and by employing 
two levels of storage with sub-iterations to converge 
the current iterate before advancing to the next time 
level 6 . 


-^-1 + #, 

At,. 


A — k + 1 

Aq,. 


( 8 ) 


The modified algorithm is conceptually represented 
by Eq. (8), identical to Eq. (4) except that sub- 
iterations on Aq k are accommodated. Equation (6) is 
modified in a similar manner. The superscript (k, 
k+1 ) signifies that the latest available iterate on 

q n + Aq k is utilized in the evaluation of r, . In 


theory, Eq. (8) is iterated on all cells until 

Aq k+1 - Aq k is converged to order At 2 . In practice 

q n+1 is updated after 10 to 20 sub-iterations on Aq . 

Algorithm Metrics: 

In most applications with crude 
initializations the point-implicit algorithm is more 
robust and efficient in reducing the residual in the 
early stages. In fact, the line-implicit algorithm will 
usually fail (increase the residuals and introduce wild 
fluctuations) when engaged too early in the relaxation 
process. However, once the solution begins to 
converge, indicated by a three order of magnitude 
decrease in the error norm, the line-implicit algorithm 
is more efficient in terms of iteration count and CPU 
time. For example. Figures 1-4 show the convergence 
history for a benchmark problem, flow over a sphere 
at V- = 5 km/s, p« = .001 kg/m 3 , T«=200K, and 
T W =500K. The convergence history is plotted as a 
function of iteration count in Fig. 1 and as a function 
of CPU time in Fig. 2 for a perfect gas model on a 
coarse grid of 30x32 starting from a uniform flow. 
The grid is not adapted. CPU times are on an SGI 
R 10000. 



Figure 1 Comparison of point- and line implicit 
iteration count for sphere benchmark case with coarse 
grid and perfect gas model. 

Subsequent grid refinements initialized from 
the previous coarse grid solution to 128 cells in the 
body normal direction gives the convergence shown 
in Fig. 3. Convergence rate slows significantly for the 
point-implicit case once the inviscid field is nearly 
converged but the viscous layer is not fully 
converged. The line-relaxation across the shock layer 
converges rapidly until round-off errors inhibit 
further error norm reduction. 

Inclusion of thermochemical nonequilibrium 
source terms in the gas model does not significantly 
change the character of convergence for the 
benchmark case on a grid of 30x64 initialized with 
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free stream conditions as shown in Figure 4. In this 
case a 7-species model for air and a fully catalytic 
wall boundary condition are specified. More frequent 
updates of the Jacobians early in the relaxation 
process are required because of the rapid evolution of 
the solution. 



Figure 2 Comparison of point- and line implicit CPU 
time for coarse grid, perfect gas, sphere benchmark. 



Figure 3 Comparison of point- and line implicit CPU 
time for fine grid, perfect gas, sphere benchmark. 



Figure 4 Typical application of point- and line- 
implicit algorithm showing CPU time for 
thermochemical nonequilibrium sphere benchmark 
case. 


Applications 

Ballutes - Spheres and Toroids: 

Large, inflatable ballutes (balloon- 
parachutes) have been proposed 7,8 as hypersonic 
decelerators for planetary aerocapture applications. 
Table 1 defines benchmark conditions for towed 
spherical and toroidal ballutes. The primary 
configuration considered here is a towed toroid, in 
which the flow disturbed by the towing spacecraft 
passes through the hole of the toroid. The toroid is 
still impacted by hypersonic flow off the tow cables, 
which are nearly in a free molecular flow regime. 
Interference effects are not included in the present 
calculations. 

Physical models in LAURA 9 have been 
updated to include chemical species expected in the 
shock layer of ballutes entering the atmospheres of 
Venus, Saturn, Titan, and Neptune - in addition to 
the models for Mars and Earth already present. The 
current model can accommodate 32 species and 57 
reaction pairs. Examples for Mars and Venus follow. 

The axisymmetric flow past a toroid at Mars 
for = 5 160 m/s and p» = 3.6 10 7 kg/m 1 (not part 
of Table 1) is presented in Fig. 5. The diameter of the 
toroid, measured from center to center of the circular 
area of revolution, is 15m. The cross-sectional 
diameter of the circle is 2.5 m. The figure shows 
vibrational-electronic temperatures in a symmetry 
plane. The computational domain includes the flow 
axis, co-incident with the axis of the toroid, the 
inflow boundary ahead of the captured bow shock, 
and the supersonic outflow boundary. The flow is 
computed with 120 cells around the toroid, and 64 
cells from the body to the outer boundaries noted 
above. The near wake is included to investigate the 
bow shock interaction along the flow axis and to 
accommodate future studies of flow interaction from 
the towing payload passing through the hole. A 
relatively high temperature interaction region is 
evident behind a Mach reflection as the intersecting 
bow shocks converge on the axis. 

The toroid for a Venus Sample Return 
mission is larger than required for the Mars Microsat 
mission noted previously. It has a more severe entry 
environment as detailed in Table 1 but the larger 
radius reduces heating rate, even with a fully catalytic 
wall boundary condition. The radiative equilibrium 
wall temperature for this case is 7 1 2 K. The peak 
shock layer vibrational-electronic temperature is 
approximately 12000 K. 
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Table 1 - Ballute Simulations 


Mission 

Venus 

Sample Return 

Mars 

Microsat 

Saturn 

Ring Observer 

Titan 

Organics Explorer 

Neptune 

Orbiter 

Atmosphere, 

Xi.« 

C0 2 = 965 
N 2 = .035 

HU 

H 2 = .963 
He = .037 

N 2 = .983 
CR,= .017 

H 2 = .80 
He = .19 
CH» = .001 

Sphere, D, m 

135. 

12. 


48. H 

100. 

V„ , km/s 

10.6 

5.54 

24.1 

8.53 

26.7 

Poo * kg/m 3 



3.1 10“ 

3.0 10 ' 

■HUHI 

C„ 

0.913 

0.99 

0.98 

0.981 ^ 

1.05 

flmax ’ W/crrr 

1.07 

1.82 

2.72 

1.93 

2.60 

T mM ,K 

677. 

774. 

855. 

785. 

845. 

Toroid, m 
Ring D / X-sec D 

150./ 30. 

15./ 3. 

120./ 30. 

52./ 13. 

100./ 25. 

V^, , km/s 

10.6 

5.49 

23.9 

8.55 

26.8 

n 

P» , kg/m 3 


MQVTjUa 

mmm 



C„ 

1.31 

1.45 

1.38 

1.39 

1.51 


1.31 

1.91 

2.85 

2.05 

2.84 

T max ,K 

712. 

783. 

865. 

796. 

864. 


from about 0.7 meters to the wall. This case was 


vlbratlaMl'Etoctranlc T*mp*r»tui» Contours computed with the line-implicit relaxation algorithm. 



z, m 


Figure 5 Vibrational-electronic temperature contours 
over toroid showing high temperature core behind 
shock-shock interaction in near wake. 

Profiles of species mole fraction across the 
shock layer are presented in Figs. 6 and 7. On this 
scale body, shock standoff distances on the order of 
several meters are computed. Significant ionization 
(20%) is also present. Note that there is no region of 
chemical equilibrium; the mole fraction of C0 2 , for 
example, does not begin significant dissociation until 
about 2 m in front of the body and shows a local 
minimum with no plateau near the boundary-layer 
edge. The catalysis model used here forces 
recombination to free stream levels at the wall; 
consequently effects of diffusion and recombination 
contribute to the rapid increase in C0 2 mole fraction 


V«nus Toroid 



Figure 6 Mole fraction of neutrals across shock layer 
of Venus toroid. 

Ballutes - Experimental Data: 

Experimental conditions 10 11 for shot 2018 in 
the T5 Hypervelocity Shock Tunnel with test gas 
C0 2 , v„ = 3.714 km/s, p„= 1.61 10 3 kg/m 3 , and 
T_= 1570K involving flow over a stainless steel 

toroid at a = 0 was simulated with LAURA. Free 
stream dissociation levels were independently 
computed 9 to be (0.32, 0.17, 0.42, 0.09) by mole 
fraction for (C0 2 , 0 2 , CO, and O) respectively. Initial 
comparisons of the experimental and computational 
results for convective heating in a C0 2 dominant 
atmosphere over a ballute are presented in Figure 8. 
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The effects of uncertainties in boundary conditions 
on predicted heating levels were not investigated. 
However, similar quality comparisons were obtained 
for shots 2013 and 2019, spanning a post shock 
Reynolds number range from 6.6 10 4 to 2.0 10. 



Figure 7 Mole fractions of ions and electrons across 
shock layer of Venus toroid. 


Run T2018 - Toroid 



Figure 8 Heat transfer over toroidal ballute measured 
in GALCIT T5 hypervelocity shock tunnel in C0 2 
and computed by LAURA. 

Thermal nonequilibrium is assumed and the 
vibrational-electronic temperature in the free stream 
is set equal to the translational-rotational temperature. 
A constant wall temperature equal to 300 K is 
specified. Boundary conditions for surface catalysis 
are varied in order to scope their influence on 
heating. In the fully catalytic boundary condition, the 
mass fraction of C0 2 is set equal to 1 at the cold 
wall. Mole fraction gradients are set to zero at the 
wall for the non-catalytic condition. In an 
intermediate approximation, it is assumed that atomic 
oxygen diffusing to the wall recombines to form 
molecular oxygen but a zero gradient is employed for 
all other species. Reasonably good agreement in 
heating is obtained with the boundary condition in 


which only oxygen catalysis is accommodated. In all 
three simulated shots, the experimental data (except 
near 0 = 90°) is bounded by the fully catalytic and 
oxygen only catalytic curves. 

Effects of Time Accurate Simulation: 

Disk and parachute-like configurations were 
also simulated in the early phases of this study but 
were not brought forward in Table 1 because of 
concerns regarding flow stability and aerodynamic 
stability. The chute configurations failed to converge 
using the constant Courant number relaxation 
scheme. Rather, they showed a periodic growth and 
collapse of the shock wave, the extremes of which 
are shown in Figure 9 for a simple Mach 6, perfect 
gas flow. Instabilities have been reported previously 
in blunt body flows with cavities 12 . While the 
magnitude of these chute instabilities was larger than 
previously reported, their occurrence was not all that 
surprising. 



Counmt nuntiMr = 10 Tima accurate 

Uiwtoady Blaady 

Figure 9 Density contours over chute at Mach 6 
showing unsteady range of motion using constant 
Courant number simulation and steady solution from 
time-accurate simulation. 

The phenomenon was interesting enough to 
revisit in order to generate a time accurate simulation 
once the capability was developed in LAURA. 
However, the time-accurate simulation shows the 
periodic motion decays to a steady state in 
approximately 8 cycles starting from t = 0. The time 
step was set at 5.23 10‘ 6 s. The reference flow time, 

D / V M , is 1.05 10" 4 s. If the constant Courant 
number simulation is initiated from a fully converged 
steady state, the steady state is retained. If the 
constant Courant number simulation is initiated from 
a “nearly-converged” time accurate solution - one 
that has seen a 4-order magnitude drop in error norm 
— the large scale periodic motion returns. Density 
contours across the shock front in Figure 9 are 
somewhat jagged because the grid had to be large 
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enough to fully contain unsteady motion and grid 
lines did not align well with the moving shock front. 

This unexpected result in which the ability 
to achieve a steady state requires temporal accuracy 
has not been observed previously in typical, convex 
blunt body geometries. However, in a recent 
hypersonic code validation study 13 , simulated flow 
over a sharp, double cone (25° / 55°) at Mach 9.5 and 
Reynolds number 264000 nf 1 in nitrogen showed 
large scale, periodic shedding of vortices in the 
separation region around the 30° compression. A 1 st 
order accurate, constant Courant number, simulation 
showed a steady flow for this same case. A 2 nd order 
simulation at a lower Reynolds number (144000 m ) 
showed a steady flow result on the same grid. 
Experimental data 14 in the Calspan - University at 
Buffalo Research Center (CUBRC) Large Energy 
National Shock (LENS) tunnel for this case (Run 24) 
from thin-film heat transfer gauges indicated that the 
flow achieved steady state in the test time. The non- 
linear minmod flux limiter in LAURA was a 
suspected source of numerical ringing within the 
recirculation region that triggered the unsteadiness. 



Figure 10 Pressure and heating over sharp, double 
cone for Run 24 of Reference 10 using time-accurate 
simulation to attain steady state. 

In light of these recent results showing that 
attainment of a steady state can depend on temporal 
accuracy in some circumstances, the sharp, double- 
cone case was revisited using time-accurate 
simulation. The time step was set at 3.65 10 7 s. The 
reference flow time, L / V M , is 7.08 10 5 s. The flow 
converged to a nearly steady state - very slight 
oscillation of streamlines persisted in the 
recirculation region, but these oscillations were 
orders of magnitude less severe than observed in the 
constant Courant number simulation. Comparisons 
with experimental data for pressure and heating are 
presented in Figure 1 0. Good agreement with the 


onset of separation is observed. Fair agreement with 
reflected wave peak pressure and heating location is 
obtained. The broad, flat pressure plateau in the 
recirculation region for 0.5<x/L<l had significant, 
varying structure in the previous, constant Courant 
number simulation. 

Both the chute and double-cone cases have a 
subsonic region contained in a concavity (rays 
perpendicular to the surface converge in a subsonic 
domain) driven by an external supersonic flow. In the 
typical grid that is highly stretched from the body a 
constant Courant number simulation means that 
physical time advances more rapidly away from the 
body than in close to the body. Perturbations 
associated with each relaxation step may tend to 
aggregate near the wall in the concavity before 
exiting a supersonic outflow boundary. Waves 
accumulate (in this non-physical temporal map), 
strengthen, and reflect back in to the subsonic domain 
in a self-sustaining but non-physical manner. This 
scenario is consistent with observed behavior in the 
two systems tested here. Certainly, not all subsonic 
domains with concavities induce this pseudo- 
unsteady behavior in constant Courant number 
simulations - there are many examples of 
compression corners with separated flow being 
computed with constant Courant number relaxation. 

It is simply noted that such flow conditions may be 
susceptible to pseudo-unsteady behavior in single 
storage level, constant Courant number relaxation 
and the overhead associated with implicit, time- 
accurate simulation may be recovered with better 
convergence. 

MSRO - Flight at Mars; 

The Mars Sample Return Orbiter is designed 
to use aeroassist to enter Mars orbit. Though Mars 
mission architectures have been revised, the MSRO 
mission investigated here called for the French 
designed orbiter to carry the US designed Earth Entry 
Vehicle (EEV) to Mars orbit. It would rendezvous 
with samples gathered in previous missions and 
launched from the surface to orbit. The samples 
would be sealed in the EEV that separates from the 
MSRO and returns to Earth. 

One of the concerns in the design is the 
requirement to protect the payload from impingement 
heating during the aeropass. Impingement may result 
from the shear layer separating from the aft edge of 
the aerobrake and turning significantly inboard. The 
shear layer turning angle has been estimated from 
previous computational and experimental studies of 
the Aeroassist Flight Experiment (AFE) vehicle from 
which MSRO design evolved. An angle-of-attack, 
a, of -4° corresponds to a 21° angle-of-attack 
relative to the base plane. At this relative angle of 
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attack the shear-layer turning angle relative to a 
backward facing normal to the base plane should lie 
between 30° and 40° (see Figure 1 1 of Reference 15). 
The payload is positioned in a region where shear 
layer impingement may be expected on the aft, 
windside corner. (Payload volume and trim angle-of- 
attack sensitivity to center-of-gravity constrain ability 
to completely avoid impingement.) Ground-based 
experimental and computational studies confirm 
shear-layer impingement on the aft corner of the 
cylindrical payload. 

Figures 11-13 show the computed, near 
wake results for the MSRO at 1 15 seconds into the 
aeropass at Mars. Flow conditions are V<* = 5223 m/s 
and p«, = 2.933 10 4 kg/m 3 at a = -4 deg. This angle 
of attack provides the strongest impingement within 
the constraints of the guidance system. The 
computational grid in the plane of symmetry is 
shown in Fig. 1 1 . Four hundred and thirty four 
simply connected, structured blocks are used in the 
near wake simulation. Inflow conditions were taken 
from a previously converged forebody solution. 



Figure 1 1 MSRO near wake mesh and surface 
pressure distribution. 

The aerobrake is removed in Figure 12 to 
reveal the recirculation region behind the cove. 
Density contours in this figure clearly show the shock 
set up over the aft corner of the payload where 
impingement occurs. Figure 13 shows a heating trace 
along the symmetry plane of the payload, 
superimposed on a global view of surface heating 
contours. Locally high impingement heating is 
consistent with observations in ground-based tests. 
Cell Reynolds number metrics at the wall are of order 
1 or less so that boundary-layer resolution should be 
adequate. However, the grid between the aft corner of 
the aerobrake and the aft corner of the cylinder 
(Figure 1 1) is not well aligned to resolve the free 


shear layer so that some numerical diffusion of the 
flow feeding the impingement is expected. 



Figure 12 Density contours and streamlines in 
symmetry plane featuring primary and secondary 
vortices behind cove of brake (removed for clarity). 



(1 000 x) / D 


Figure 13 Trace of heating on symmetry plane of 
payload showing maximum value at impingement 
and local maximum under secondary vortex in cove 
region. Heating contours are on aftside of aerobrake. 

MSRO - CF^tunnel; 

Wind tunnel tests of the MSRO 
configuration were executed in the Langley CF 4 
tunnel in order to study impingement heating levels 
in a gas with a relatively low ratio of specific heat - 
as will be encountered in flight 16 . The configuration 
is a scale model of the flight vehicle simulated above 
but a new grid (Figure 14) was generated for this 
study to simplify blocking topology (70 simply 
connected blocks) and enable automatic grid 
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adaptation into the far field. The new grid has finer 
resolution in the free shear layer but is still not well 
aligned with streamlines approaching the aft corner. 



Figure 14 Grid over MSRO for CF4 wind tunnel 
study with surface pressure distribution. 

Simulated flow conditions are: = 907.4 

m/s, p„= 1.53 1 O' 3 kg/m 3 , and T M =201.1K. Fully 
laminar flow is expected at this condition based on 
observed trends in the experiment as a function of 
Reynolds number 1 ^. Heating along the leading edge 
of the cylindrical payload (Figure 15) is non- 
dimensionalized by a reference stagnation point 
heating level defined by the method of Fay and 
Ridell. This reference heating is about 10% larger 
then measured stagnation point heating for a related 
case - only a few forebody measurements were made 
because the study was focused on payload heating 
issues. Comparisons between experiment 
(thermophosphor technique) and computation were 
generally within 8% on the forebody in this related 
case. The abscissa of Figure 15 refers to the running 
length along the cylinder non-dimensionalized by 
aerobrake reference diameter, originating at the aft 
corner and running forward towards the aerobrake. 

The experimental data (symbols) indicate a 
glancing impingement, in which heating is highest at 
the corner. The LAURA results indicate a stronger 
impingement, in which a plateau level of heating is 
predicted, and a subsequent heating spike at the 
corner (x/D = 0) where there is a rapid thinning of the 
boundary layer as flow expands into the wake. A grid 
refinement (factor 2 in the streamwise and 
circumferential directions) provides better resolution 
of the shear layer but shows essentially no change 
(less than 3%) in surface heating. A heating spike 
near x/D=0.3 is grid induced - 4 blocks fan out from 
this location to define the cove region and the 
“streamwise" coordinate of the two central blocks in 


the fan spans the boundary layer with inadequate 
resolution. 



Figure 15 Impingement zone heating on cylinder 
along plane of symmetry. 

The flight simulation of heating from Figure 13 
is rescaled and plotted in Figure 15. In this case, q re f 
is the computed heating level at the stagnation point 
on the forebody. At an earlier point in the trajectory 
(t=70s - not plotted here) the same trend is observed 
but it scales lower such that the maximum value is 
q/q ref =0.3. The grid along the cylinder was coarser in 
the flight case but did not have a singularity. There is 
no plateau in heating - possibly because the grid is 
too coarse to support one. Still, the general trend is 
that LAURA predicts a larger area of high heating 
surrounding impingement on the given grids than 
measured in CF4. Additional simulations are planned 
with re-aligned grids and, possibly, application of slip 
boundary conditions in attempt to resolve differences 
in the CF 4 case. 

Concluding Remarks 

Point- and line-implicit relaxation 
algorithms within Program LAURA for aeroassist 
applications are defined. The strategy of opening 
with a more robust point-implicit solver and closing 
with a more aggressive and efficient line relaxation 
was applied to every problem in this paper. 
Convergence metrics are documented for a model 
problem. Implicit wall boundary conditions for the 
fully coupled, thermochemical nonequilibrium 
system are required to make the line relaxation 
strategy work. These conditions are defined herein. 

Hypersonic flows with subsonic regions 
bounded by a local concavity were observed to 
support periodic, unsteady behavior using a constant 
Courant number relaxation but recovered a steady 
flow behavior using time-accurate simulation. The 
behavior is thought to emerge from the temporal 
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implicit method abruptly diverges just as a fully 
converged state appears to be at hand. This behavior 
is observed even when the implicit treatment of mole 
fraction in the ghost cell uses the fully catalytic 
approximation - the demise just takes longer. Such 
behavior is indicative of a singularity in the 
coefficient matrix associated with the converged, 
non-catalytic state but this conjecture has not yet 
been confirmed analytically. 
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